Screening of a hypercritical charge in graphene 
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Screening of a large external charge in graphene is studied. The charge is assumed to be displaced 
away or smeared over a finite region of the graphene plane. The initial decay of the screened 
potential with distance is shown to follow the 3/2 power. It gradually changes to the Coulomb law 
outside of a hypercritical core whose radius is proportional to the external charge. 
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Recent discovery of graphene - a two-dimensional (2D) 
form of carbon [l| - brought an exciting link between 
solid-state physics and quantum electrodynamics (QED). 
The half-filled 7r-band of graphene has a relativistic mass- 
less Dirac spectrum e = ±?iu|k| where e > for the 
electrons and e < for holes, k is the deviation of the 
quasi-momentum from the Brillouin zone corner, and 
V « lO^m/s. The role of the fine-structure constant is 
played by the dimcnsionlcss parameter 



e^/Khv, e'^/hv fii2.2, 



(1) 



where k is the dielectric constant at the interface of sub- 
strate and vacuum. For conventional Si02 substrates 
K « 2.4; hence, Coulomb interaction is strong, a ~ 1. 

In this work we consider the problem of screening of 
a Coulomb potential Vq = eZ/ nr that can be induced in 
graphene by a group of charged impurities in the sub- 
strate, by a nearby gate, or by a cluster of dopants. 
This problem is important for a number of properties of 
graphene nanostructures, including transport 
local gating B, B, 0|, controlled doping dOjIlll) and 
chemical sensing ll a. Not surprisingly, it has attracted 
much attention [all, i, E [liJllKa^ [ll [H. In 
particular, it has been noted 

mm \lm that at half- 
filling the problem of a Coulomb charge in graphene has 
an interesting parallel with that of a hypothetical su- 
percritical atom with Z > Hc/e^ « 137. For such an 
atom, the standard solution [20| of the Dirac equation 
breaks down and a physically acceptable atomic struc- 
ture is obtained only after accounting for a finite radius 
of the nucleus 2l| . This structure is characterized by a 



vacuum reconstruction: a certain number of electrons is 
spontaneously created (liberating positrons), they bind 
to the nucleus, and render it subcritical. In graphene the 
critical charge l^, T^iEillSl -^c — l/2a is much smaller 
than in QED; hence, solid-state analogs of supercritical 
atoms may be realizable even at Z ^ 1. 

According to all prior investigations, screening proper- 
ties of an undoped graphene resemble those of a dielec- 
tric: the screened potential ^ of a supercitical charge has 
been argued not to deviate much from the Coulomb law, 

e F{r) 



Here F{r) is a slow logarithmic function. Such a conclu- 
sion follows from the standard linear response theory — 
Random Phase Approximation (RPA) [J — and was 
supposedly confirmed by calculations within the Thomas- 
Fermi (TF) method that is able to go beyond 
the linear response. Below we re-examine these conclu- 
sions for the case of a hypercritical charge Z 3> 1, which 
lets itself to a controlled treatment and adds new physics. 
Without loss of generality we assume that the external 
charge attracts electrons, Z > 0. 

Since it was not always made clear previously, we em- 
phasize that the problem is ill-defined unless one explic- 
itly regularizes the strong Coulomb singularity at the ori- 
gin. This is as crucial as introducing a finite size of a 
nucleus in QED. Therefore, the charge Z must be ei- 
ther displaced away from graphene plane by some dis- 
tance d or spread over the area of some radius in this 
plane. In order to deal exclusively with Dirac fermions 
the smearing parameter maxjd, ro} must exceed aa.\fZ 
where a = 2.5 A is the graphene lattice constant; other- 
wise, the quasiparticle energy shift due to the potential V 
would exceed the modest energy separation 4eV e^/a 
of the Dirac point and the nearest cr-bands [l^. These 
other bands would then also need to be included, leading 
one to a three-dimensional (3D) problem that has little 
to do with special properties of graphene. 

Our main result is that the induced 2D electron density 
and the screened potential have the form 



71 (r) 



1 n 



V{r) 



ri = 2a^Zd (3) 



y(r) = 



Kr 2a^ 



(2) 



in the range of distances maxjd, ro} ^ r <C ri. Based on 
electrostatics, the law ([3]) is robust and universal. How 
does then one reconcile it with Eq. ([2])? As we clarify 
below, the situation is as follows. In the strongly inter- 
acting case, a ~ 1, Eq. ^ controls the entire supercrit- 
ical core, i.e., the circle around the origin where the net 
charge exceeds Zc- This fact has eluded previous studies. 
However, if a is small, the domain of validity of Eq. ([3]) 
narrows down, opening up a window where Eq. ([2]) is 
realized. Although current experiments are not in this 
regime, small a can be achieved using large k substrates, 
e.g., Hf02 3 or simply liquid water, k ~ 80. 
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The three-line derivation of Eq. ([3]) can be given if, as 
discussed above, the charge Z is point-hke but removed 
from the graphene plane [l3| by an appropriate distance 
d. The key idea is that if we treat the graphene sheet as 
a perfect metal, then classical electrostatics dictates that 
the induced charge density is given by n = ?ici, where 



nci(r) 



1 



Zd 



1 



27r (r2 + d2)3/2 47rQ-2 (^2 ^2)3/2 



(4) 



At r ^ d we get the first formula in Eq. ([3]). To derive 
V{r) we employ the TF approximation, 



fi[n{r)] ~ eV{r) = 0. 



(5) 



Combined with the formula for the chemical potential, 

^(n) = sign{n)^/Trhv\n\^^^ , (6) 

specific for the 2D Dirac spectrum, it yields the second 
formula in Eq. concluding the derivation. The rest 
of our paper is needed mainly to explain why the above 
reasoning is correct, why Eq. ([3|) is completely general 
rather than restricted to the case of a remote charge, and 
finally, where the room may still exist for the differing 
predictions for n and V advocated in Refs. 

First, let us clarify why it was legitimate to approxi- 
mate the density response of graphene — a complicated 
quantum system — simply by that of an ideal metal. 
The reason is this. At r <C ri the local screening length 
rs — (k/ 27re2)(d/i/dn) ~ Q;~-'^|n|~^/2 jg niuch smaller 
than the characteristic scale max{r, d} over which the 
potential V{r), or equivalently, the effective background 
2D charge density nci(r) vary. Therefore, the unscreened 
charge density, cr{r) = rici(r) — n{r), is smaller than the 
background one, nci(r), by some large factor related to 
the ratio of rs and max{r, d}. [The precise relation is 
expressed by Eqs. (fTTj) and (fT2|) below.] 

The next step is to explain why or rather where the 
TF approximation can be trusted. This is determined 
by the conditions that max{r, d} exceeds the local Fermi 
wavelength Xpir) ^ n^^^^{r). For a ~ 1 we can use 
n(r) from Eq. ([3|) to write this condition as r < r2 = Zd. 
Thus, for a ~ 1 the domains of validity of the TF and the 
perfect screening approximations coincide, ri ^ r2. At 
r <^r2 all corrections to Eq. both smoothly varying 
with r and Friedel oscillations [16jjl24| are subleading. 

Let us briefly discuss the nature of screening at r > r2 
where the TF approximation breaks down. Define Q{r) 
to be the net effective charge inside the circle of radius r. 



Q(r) = / 2na{r'ydr' . 



(7) 



At r — r2, Q drops to a number of the order of the 
critical one Zc ~ l/2a. Consideration of screening 
now requires a detailed analysis of the eigenstates of the 



Dirac equation [l^ 12 j in the potential created by 
the charge (5(r2). According to Ref. some amount 
of charge, in fact, exactly the critical one remains un- 
screened: Q{oo) = Zc- The saturation of Q at this value 
occurs near a certain r = r». However for a ~ 1, 
and r2 must coincide up to a factor of the order of unity. 
Thus, Eq. ([3]) governs the entire supercritical core except 
perhaps a non-parametrically wide outer region r ~ r2 
where a more complicated dependence [l6| may apply. 
At even larger distances the potential V{r) follows the 
RPA prediction 



V{r) ~ eZc / er , r » r2 , 



(8) 



where e = k[1 + (7r/2)a] is the RPA dielectric con- 
stant 0, • A more careful examination of the behav- 
ior of V(r) at such r requires accounting for the infrared 
renormalization of a (which enters e) |l4l . lol . 27 1, that is 
not directly related to the problem at hand. 

Let us return to the analysis of the supercritical region 
and show how to refine our results by computing correc- 
tions to Eq. For this we complete the set of the TF 
Eqs. ([5]) and ^ by adding another one for V{r): 



K fd'^r'air') f , , , , f 
-e^^JT~^^ y dgJo(gr)a(<;) = J 



g{s)ds 



(9) 



where Jo(^) is the Bessel function [25| and cr, aptly 
parametrized by a'(q) = ds g(s) cos qs j26| , is the 2D 
Fourier transform of a. Inverting the last equation of ([9]), 
we get g(u) = {2K/TTe)(d/du) J^^ V{s)sds/\/ v? — and 



Q{r) = Q((X3) 



2 K 

TT e^ 



idu 



r2 du 



eV{s)sds 



(10) 



The leading correction to the perfect screening can 
be obtained by substituting /Lt[nci(s)] in lieu of eV(s)^ 
cf. Eq. (O. The resultant expression is cumbersome, and 
so we quote only the limiting forms: 



nci{r) 




r <^d, 



(11) 



r4(5/4)./-, d«r«ri, (12) 



where the Gamma-function [25|] r(5/4) « 0.906. In 
agreement with the above physical argument, the devia- 
tion from the perfect screening at all r <C ri is small. 

These analytical predictions were verified by numerical 
simulations. To this end we solved the TF equations 
([9]) inside a finite square of the 2D plane. The integrals 
were replaced by discrete sums over a uniform 256 x 256 
grid defined therein and the periodic boundary conditions 
were imposed. The solution for n(r) and V{r) was found 
by a standard iterative method, using underrelaxation to 
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ensure convergence. As shown in the inset of Fig. [TJ the 
analytical and the numerical results agree extremely well 
for a suitably large hypercritical charge o?Z = 20. 

Small-a regimes. — Let us now show that Eq. Q can 
be reconciled with our theory under the condition a <^ 1, 
i.e., k:s> 1. In this case there is a gap between the above 
defined characteristic lengthscales ri and r2. This gap is 
filled by an additional regime where the TF approxima- 
tion regime is still valid but the screening is ineffective. 

To see that consider first moderately small a, such that 
<C a <C 1. Since the screening is weak, Eq. (jTU]) 
is no longer convenient. Instead, the derivation of V 
and n can be done along the lines of Ref. but with 
several important refinements. First, we trade the two 
last equations of ^ for 



-V{r) 



^Kr-^][n^,{r')-n{r')], (13) 



where K{z) is the complete elliptic integral of the first 
kind [2^. Next, we treat Eq. ^ as the definition of yet 
unknown function F and use Eqs. ([5]) and ^ to obtain 



n{r) = F^{r) / AiTa^r'^ 



(14) 



Taking the limit d ^ at fixed ri, we get the equation 

oo 

F{t)= [ du[eit~u) + <j){u-t)][e-'' - F^{u)], (15) 



where t — ln(r/ri) and 0(t) is the unit step- function. 
Function <j){t), defined by Eq. (12) of Ref. [11], has the 
following properties: it is a logarithmically divergent at 
t = 0, is exponentially small at \t\ ^ 1, and satisfies 
J (t){t)dt = In 4. It is easy to see that at large negative 
t, we must have perfect screening, F'^{t) ~ e~*. The 
asymptotic behavior of F at large t > can be deduced 
by replacing (t){t) with (ln4)(5(t) 28]. After this, we can 
differentiate the integral equation (fT5| to get 

F"^ - (21n4)lni^ = ln(r/ri) +c, c = const , (16) 

where we returned to the original linear coordinate r. 
The direct numerical solution of Eq. (|15p shows in excel- 
lent agreement with Eqs. ^ and Eq. if the constant 
c is set to 0.6, see Fig. [TJ At r ~ ri, this solution crosses 
over to the strong screening regime, Eq. p2p . 

The range of r ^ r2 where Eq. (jl6|) is valid is again 
determined by the condition r 3> n^^/^(r), which yields 



r2 ~ ri exp(l/ 2-/?? 



a 



(17) 



At r ^ the screened potential is given by Eq. 

Consider even smaller a, such that 1/Z <C a ^ 1/^/Z. 
Here the Coulomb interactions are so weak that the 
smearing of the external charge is no longer necessary: 
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FIG. 1: (Color online) Main panel: Density profile in the 
limit of d — > at fixed ri. The quantity plotted on the ver- 
tical axis is 47ra^rin(r) — F'^(r)rf/r^ . The thin black line is 
the numerical solution of Eq. (|15p : the red dashed line is the 
perfect screening, F — ri/r; the thick cyan line is for F from 
Eq. (|16p with c = 0.6. Inset: expanded view of n(r) inside 
the hypercritical core. The thin black line and the red dashed 
line have the same meaning as before; the dots correspond to 
an analytical formula whose limits are given by Eqs. (|ll|l and 
CI. 



the "dangerous" region r < aa\j Z is smaller than the 
lattice constant. In addition, the domain of the perfect 
screening, r < ri , which is the region of validity of Eq. ^ 
disappears. (Weak interactions entail poor screening.) In 
this case c — * {2a'^Z)~^ -\- ln(ri/a), so that the solution 



V{r) 



eZ 



1 



(18) 



advocated in Ref. [15| actually applies, at ln(r/a) <l/a. 

In-plane charge. — In the concluding part of the pa- 
per we wish to return to the structure of the hypercrit- 
ical core and to show that Eq. ^ remains valid if the 
charge Z resides within the 2D plane. To gain some intu- 
ition consider first an artificial scenario where the exter- 
nal charge is highly localized yet the tr-bands of graphene 
can be neglected. In this case the maximum possible elec- 
tron density (relative to that of the half-filled 7r-band) 
is ^^lnax = 2/-\/3a^. This density is indeed reached at 
r smaller than some radius 6 as a result of attraction of 
electrons to the hypercritical charge Z. At r > 6, electron 
density is gradually decreases, which can be thought of 
appearance of "holes" at the top of the conduction band. 

Incidentally, the charge profile of these holes within the 
perfect screening approximation is known exactly. It can 
be read off the results of Ref. [2§] where the structure of 
a depletion region around a disk of a negative charge in a 
semiconductor was studied. For a high density of the ex- 
ternal charge these authors found that h = (Z/27rnoo)^^^, 
where rioo is the uniform electron density far away form 
the depletion. They also found [1^ that the density pro- 
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file at large r is given by n{r) = Uoo — {Zb/2Trr^) atr^b. 
Adopting these results to our problem, we get 

n{r) = Zb/2nr^ , r ■> b = (Z/2™„,ax)'/' , (19) 

leading to Eq. © with ri = 2a^Zb - aZ^/^ for a - 1. 

Consider now a more realistic setup where the external 
charge ricxt (r) is distributed over a disk of radius rp ^ 
a\fZ . Then n < n(0) ~ Z/tttq rimax, so that cr-bands 
bands can indeed be disregarded. Let us show that 



n(r) 



ZTs 



ro ^ r <C ri = 2a Zrg 



(20) 



where the screening length Vs ^ l/ay/n{0). 

Based on the near-perfect screening framework used 
in the first part of the paper [and justified a posteriori 
by Eq. ([20]) ] we can claim that V{r) is substantial only 
in the region r < tq and is greatly reduced at r > tq. 
This implies that the Fourier transform of V is nearly 
wavevector-independent over a range of q, 

Viq) = CiZers/fi + 0{Z^/^) , r^^ « q « r^^ . 

The first term, with ci ~ 1, follows from Eq. ([5]). 
In turn, the Fourier transform of the charge density, 
a{q) = hcxtiq) — n(q) = V{q)/{2Tre/Kq), that produces 
this potential behaves as a (g) ==C2Zr,g+C'(Zi/2), where 
C2 ^ 1. After the inverse Fourier transform, the net 
charge density (7{r) is seen to be dominated by the term 
—C2Zrs/2Trr^ at <C r <C ri. Since ricxtif) — for such 
r, this term is entirely due to n, proving our statement. 

In summary, we considered the problem of nonlinear 
screening of a large charge by the massless electrons in 
graphene. The consistent formulation of the problem re- 
quires the charge to be either displaced from the graphene 
plane or to be spread over a disk of finite radius rp . In 
both cases the screening is nonlinear within a region of 
a parametrically large radius ri. In the interval between 
ro and ri the screened potential decays as 1/r^^^. Our 
results are relevant for current and future experiments 
that involve local charging or doping of graphene. Thus, 
if small a can be achieved experimentally, it may be pos- 
sible to verify the predicted crossover from Eq. ([3]) to ([2]) 
and finally to (fT8|) by using scanned probe techniques [9| . 
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